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SUMMARY 


This report describes the second phase of an analytical ef- 
fort to model the effects of low velocity transverse normal im- 
pact on laminated composite plates. The methodology utilized 
consisted of a transient dynamics finite element analysis with 
capabilities for composite material characterization and stress 
predictions. A shear flexible plate bending element with mate- 
rial bending/extensional coupling capabilities was developed for 
this effort. The contact effects between the impacting mass and 
composite plate were modelled utilizing a nonlinear contact spring 
similar to Hertzian contact. Provisions were included for incor- 
porating the effects of both ply damage and delaminations through 
modifications to the element stiffness matrices. 

An impact analysis was performed to simulate the impact of 
a steel sphere on an eight ply, quasi-isotropic graphite /epoxy 
plate. The response of the plate and impactor were modelled 
through 100 ysec of the impact event. Damage predicted included 
both matrix cracking and fiber breakage within the composite plate. 
Damage was shown to be a function of high flexural modes present 
early in the impact event. 
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INTRODUCTION 


With the ever increasing use of laminated composites in struc- 
tural applications, an interesting phenomenon has become apparent. 
This phenomenon concerns the real possibility of invisible damage 
within a composite structure caused by low velocity, low mass im- 
pacts. This type of loading environment is most easily envisioned 
as the impact of a dropped workman's tool on a structure or, per- 
haps, runway debris ejected onto a structure. 

The primary concern related to this type of loading environ- 
ment is the introduction of performance degrading damage within 
the composite structure. This damage may not be visually apparent 
during subsequent inspections of the structure and, hence, the per- 
formance degradation will also not be immediately apparent. 

The subject of impact related phenomena has been studied by 
many investigators utilizing many different approaches. Much of 
this work has been related to ballistic type impact and, hence, is 
not applicable here. In ballistic impact, the velocities involved 
are high enough to promote large stress wave propagation effects. 

Prior work in the analysis of the low speed impact problem 
has established that it is reasonable to neglect the stress wave 
propagation problem and to focus on the transient structural dynamic 
approach. Different approaches have appeared in the literature to 
combine contact effects with dynamic effects. The Hertz contact 
problem has been extended to the problem of dynamic contact and 
also to the problem of contact of anisotropic bodies (see ref. 1). 
These approaches treat impact with a semi-infinite target. In the 
present case, one is concerned with a target in which the dynamic 
response of the target is important in the sense of transient 
structural motion rather than material displacement. This problem 
appears to have been addressed first by Timoshenko (ref. 2) , as de- 
scribed by Goldsmith (ref. 3) . Timoshenko studied the problem of 
the impact of a beam where the contact between the bodies was gov- 
erned by Hertz's law for contact deformations. Karas (ref. 4) ex- 
tended the Timoshenko approach to the study of plate impact (see 


ref. 5) . Moon (refs. 6 and 7) has utilized the Hertz impact theory 
in combination with a Mindlin plate theory (ref. 8) to model a sim- 
ilar approach for impact of plate structures. This particular ap- 
proach yields a nonlinear mathematical problem and extended numeri- 
cal analysis is required to obtain solutions. The procedure is suf- 
ficiently complex to motivate consideration of alternate approaches. 

Two such approaches are based on simplifications of the con- 
tact force analysis. In one case, it is considered that the impact 
takes place during a time period which is very short compared to 
the period of the first natural frequency. In this case, it appears 
reasonable to regard the impact as having imparted an impulse lo- 
cally to the plate and to then study the dynamic response of the 
composite plate target to that impulsive loading. This approach 
(see ref. 9) is appropriate as the structural stiffness increases 
and the impacting mass decreases. 

Another line of approach initiated by Clebsch (ref. 10) as de- 
scribed by Goldsmith (ref. 3) assumes that upon impact, the projec- 
tile moves with the plate and that the velocity of the projectile 
becomes an initial velocity condition. Thus, the analysis is the 
structural dynamic response of the plate with the attached mass. 
McQuillen et al. (ref. 11) applied this approach to a beam. They 
minimized some of the numerical problems by considering the con- 
tact zone between projectile and target to have finite width. This 
approach tended to minimize the contribution of the higher frequency 
modes and, thus, numerical procedures were more successful. How- 
ever, even with these assumptions, the work of reference 11 shows 
that modes of vibration other than the fundamental mode can be ex- 
cited by impact, particularly if the striker mass is small. This 
approach, which is expanded somewhat in references 12 and 13, is 
being utilized in the current effort. 

In the present study, the primary emphasis has focused on the 
prediction of damage initiation and propagation during the impact 
event and subsequent dynamic plate response. The difficulty posed 
here is that any induced damage will alter the plate stiffness lo- 
cally and, hence, affect the subsequent dynamic response. Thus, 
the closed form analytical approach taken in references 11-13 is 
not applicable. 
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In the first phase of the current effort (ref. 14) , a classical 
thin shell finite element transient dynamics analysis (ref. 15) was 
coupled with lamination theory to predict the flexural response of 
composite plates subjected to transverse impact. The impact was 
assiamed to be perfectly plastic with the impactor mass lumped at 
specified nodes within the finite element model. Both layer and in- 
terlaminar damage were predicted using a stress based failure cri- 
terion. The effects of impact induced damage were included by mod- 
ifying the local, elemental stiffness matrices to reflect the loss 
of ply stiffness. No stiffness reductions were possible for inter- 
laminar damage because the model did not include the effects of 
shear deformation. 

In this second phase, the transient structural dynamics model 
with stiffness modifications to account for damage has been retained. 
The model has been extended to incorporate nonlinear contact effects 
between the impactor mass and the plate structure. This effort re- 
flects the desire to model the energy transfer mechanisms and energy 
loss which occur at the point of impact and utilized the work of Sun 
and Yang (refs. 16 and 17) . 

In addition to the contact effects, the structural dynamics 
model has been extended to include the effects of shear deformation. 
This required a new shear flexible plate bending finite element de- 
rived from the work of Hinton et. al. (ref. 18) and transverse 
shear constitutive relations developed by Cohen (ref. 19) . 

The computerized procedure has been utilized to predict the 
response of a laminated plate subjected to a low velocity impact 
event. Structural response including damage initiation and growth 
has been predicted and is documented within this report. 
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ANALYTICAL METHODOLOGY 


The transient dynamic analysis of laminated composite plates 
subjected to low velocity impact can be broken into three inter- 
related procedures: stiffness formulation, displacement solution 

and laminate stress/failure analysis. The stiffness of the model 
depends on the failure status of the laminate; the displacement 
solution depends upon the stiffness of the model; and, the state 
of stress in the model depends upon the displacement field. The 
stress state defines the failure status of the laminate which re- 
defines the stiffness of the model, and the process continues. 

The analysis begins with the formulation of the stiffness 
for an undamaged laminate, using conventional finite element tech- 
niques. The finite element, which has been developed for use in 
CLIP II, is an eight noded isoparametric quadrilateral with cubic 
displacement functions. This element possesses both plate bending 
and membrane displacements, and includes shear deformation effects 
and the ability to handle materials with bending-extensional cou- 
pling. These features make it an ideal element for the study of 
thin and thick laminated shells which need not be balanced or sym- 
metric. The theoretical details of the formulation of this ele- 
ment are included in Appendix A. 

Being a dynamic analysis, the mass and damping matrices for 
the systems are also required. A diagonal (lumped) mass matrix 
is used, which includes rotary inertia terms. As the analysis 
proceeds, the stiffness of any individual element may change to 
reflect the effects of ply damage and delaminations; the mass 
matrix, however, remains unchanged throughout. It is assumed that 
the damping matrix is a linear combination of the mass and stiff- 
ness matrices (Rayleigh damping) and will therefore change in ac- 
cordance with changes in stiffness. 

To completely describe the mass and stiffness of the dynamic 
system, the mass and contact stiffness of the impactor must also 
be known. The impactor is treated as a liamped mass that is tied 
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to the plate at a single node by a nonlinear contact spring. • The 
nonlinear contact spring contains hysterisis, and its stiffness 
vanishes if the plate and impactor have separated. Details of the 
contact spring are given in Appendix B. 

With the mass, damping and stiffness of the system known, the 
equations of motion may be integrated to yield nodal displacements, 
velocities and accelerations as functions of time. CLIP II per- 
forms this integration momerically , using the fourth order Runge- 
Kutta-Gill method. The R-K-G method is a very accurate algorithm, 
but has the disadvantage of being sensitive to all dynamic modes 
of the system. This may require that a relatively small time step 
size be used to insure stability. 

The procedure is started by specifying the initial conditions, 
both displacements and velocities, of the system. CLIP II has pro- 
visions for specifying nonzero initial conditions for the plate and 
the impactor. Note that gravity is not included in the model, 
therefore an impactor that is given an initial displacement with no 
initial velocity will not "fall" on the plate. 

A restart capability has also been included in the code. When 
a restart analysis is performed, the last solution is read from 
storage and used as the initial conditions for subsequent solutions. 
The time increment used in the restart analysis may be different 
than that used in the original analysis. This capability allows 
the analysis to be stopped and restarted with a larger time incre- 
ment after higher frequency modes have damped out of the problem. 

An option to correct the dynamic solution after a specified 
number of time steps has been included. This is done using Ham- 
ming's fourth order corrector algorithm. This algorithm is an im- 
plicit scheme which uses information from the last three time steps 
to correct the solution at the present time. Being an implicit 
scheme, it must be used iteratively, which may become very time 
consuming. It is typically more cost effective to improve the ac- 
curacy of the results by reducing the size of the time increment 
in the R-K-G procedure, rather than evoking this correcting algo- 
rithm. Hamming's method also requires that the time increment be 
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constant for the three time steps prior to the current time step, 
which may prohibit its use on some solutions. The R-K-G method, 
Hamming's algorithm, and their implementation in CLIP II are dis- 
cussed in further detail in Appendix C. 

The displacement solution at any time step may be used to per- 
form a laminate stress analysis on an element by element basis. 

The strains at the centroid of an element are found by averaging 
the strains calculated at the Gauss points for the element in ques- 
tion. These are computationally efficient points to use, since the 
strain-displacement relationships have already been evaluated there 
during the stiffness formulation. These strains and curvatures are 
found in the global coordinate system, and must be transformed into 
material coordinates on a ply by ply basis. 

Stresses are calculated at the bottom, mid-plane and top sur- 
faces of each ply. These stresses are then used to predict ply 
damage or delaminations between plies. The failure criteria used 
in CLIP II are listed in Appendix D. These criteria predict both 
the presence of failure and the mode of failure. 

The mode of failure determines how the element stiffness matrix 
will be modified to reflect the effects of the predicted damage. 

Ply failures and delaminations are handled in radically different 
ways. Ply failures are incorporated by modifying the material prop- 
erties of individual lamina, whereas delaminations are included by 
modifying the fundamental strain-displacement relationships for the 
element. 

The way in which the lamina properties are changed depends up- 
on the type of ply damage. If the ply damage includes fiber fail- 
ure, then the entire ply is removed from the laminate. If matrix 
failure is the only failure mode, it is assumed that the fiber can 
still carry load; thus, only the transverse properties of the ply 
are removed. These laminate modifications apply only to the spe- 
cific element under consideration. Damage in one element will, 
therefore, not affect the stiffness of other elements. 
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Delaminations pose a more difficult problem. CLIP II incor- 
porates delaminations by modifying the strain-displacement rela- 
tionship to reflect the effects of singular forces in the plane 
of the delamination. The magnitude of this force can be expressed 
in terms of the displacement functions and the constitutive rela- 
tions for the laminate, which allows it to be incorporated in the 
finite element scheme. The details of this approximate procedure 
are included in Appendix E. Unfortunately, the formulation used 
for incorporating delaminations in the analysis places some severe 
restrictions on the geometry of the element. Only rectangular 
elements may be used, and they must be oriented as shown in fig- 
ure 1. 

Upon completion of the laminate analysis, the stiffness of 
the model is modified to reflect any damage that may have been 
predicted. The dynamic analysis then proceeds at this reduced 
stiffness until further damage is predicted by the laminate analy- 
sis and additional stiffness reduction is performed. This proce- 
dure continues until the desired number of time steps has been 
reached, or an element suffers complete failure. 
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IMPACT ANALYSES 


The impact analysis code described earlier has been utilized 

to predict the response of 10.16 cm x 15.24 cm rectangular clamped- 

clamped graphite epoxy plate. The laminate is comprised of T-300/ 

5208 material arranged in a [45/0/-45/90] layup. The material 

s 

elastic and strength properties utilized are listed in table 1. 

The laminate total mass was 25.3 g with a total thickness of 
0.1057 cm. 

The impact event modelled corresponds to a 1.59 cm steel 
sphere (16.45 g) impacting the plate centrally with a velocity of 
9.4 m/s. The various coefficients required to define the response 
of the contact between the plate and impacting sphere were deter- 
mined utilizing data from references 16 and 17 and are listed in 
table 2. 

CHARACTERIZATION OF IMPACT RESPONSE 


An initial analysis was performed without stress calcula- 
tion in order to determine both the characteristics of the impact 
event and verify the workings of the CLIP II code. This analysis 
also provided information relating to the time step required for 
numerical stability of the R-K-G solution procedure. 

The finite element model used in the initial analysis is 
shown in figure 2. This model is rather crude but still provides 
a reasonable displacement response and yielded the required in- 
formation. Since no stress analysis was included, the analysis 
yielded the elastic response of an undamaged plate. 

The dynamic response was initiated with an integration time 
step of 0.1 Msec. After approximately 10 ysec, the time step was 
gradually increased until a maximum time step of 0.25 ysec was 
found. Above 0.25 ysec, the solution procedure became numerically 
unstable and divergent. This very short time step requirement 
stems from the sensitivity of the solution procedure to all possi- 
bly dynamic modes present in the model. In the solution where the 
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integration time step was larger than 0.25 ysec, the rotational 
degrees of freedom appeared to become unstable first indicating 
either high bending modes and/or transverse shear effects were 
causing the difficulty. When the time step was specificed to be 
0.25 ysec, no difficulties were encountered and the plate response 
was monitored out to 750 ysec. This period of time was long 
enough to include the maximum lateral deflection of the plate 
and hence characterize the plate response to the impact event. 

The lateral deflection of the node at which the impact occurs, 
and the displacement of the impactor have been plotted as functions 
of time in figure 3. Notice that the impactor follows a smooth, 
almost sinewave like curve, while the center of the plate travels 
an irregular path. This implies that higher dynamic flexural 
modes are combining with the fundamental flexural mode to define 
the response of the plate. 

Although the mass of the entire plate is of the same order 
of magnitude as the impactor, the effective mass of the portion of 
the plate which undergoes substantial motion is far less than that 
of the impactor. Thus, by conservation of momentum, it is possi- 
ble for the plate to achieve a velocity which will allow it to 
separate from the impactor. 

Separation of the plate and impactor occurs three times dur- 
ing the time interval analyzed. In figure 4, separation of the 
plate and impactor is seen to occur from approximately 100 to 170 
ysecs, again from approximately 280 to 440 ysecs and finally from 
540 to 560 ysecs. During each of these intervals, the plate moves 
under the action of inertial forces alone since the contact force 
between the plate and impactor has vanished. 

Referring back to figure 3, the plate/ impactor separation 
during the period 280 to 440 ysecs is readily apparent. During 
the other two separation intervals, the total separation distance 
is smaller than the permanent deformation which has occurred in 
the plate due to the contact of the impactor. Hence, the loss 
of contact cannot be seen directly. 
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When contact is lost, the impactor is in "free fall" (con- 
stant velocity) , until contact is regained, once again causing 
the impactor to lose velocity. The impactor continues to slow 
down until approximately 650 ysec into the event. At this time, 
the plate begins to act like a compressed spring which releases 
its energy by pushing the impactor back in the direction from 
which it came. The velocity of the impactor passes through zero 
and changes sign. It will continue to pick up speed until it is 
thrown free from the plate, ending the event. 

The displacement response very early in the impact event is 
depicted in figure 5. Displacements through the center of the 
plate along the longer dimension are shown at 25 ysec to 100 ysec. 
Throughout this period of time, the plate response is easily seen 
to be dominated by flexural modes other than the fundamental mode 
of the clamped-clamped plate. At 100 ysec, the maximum plate de- 
flection is approximately 75% of the plate thickness. This would 
indicate that the flexural response is still the dominant load 
carrying mechanism at 100 ysec. Referring to figure 3, the max- 
imum displacement is seen to approach the plate thickness at ap- 
proximately 125 ysec. As the maximum displacement becomes larger 
than the plate thickness, the stiffness of the plate will increase 
as large deflection membrane action becomes significant. 

Since the analysis performed is strictly linear, these ef- 
fects are not included. Hence, beyond approximately 125 ysec 
the accuracy of the predicted response must be considered suspect. 
The trends which have been shown are believed to be representa- 
tive of the impact event. The multimode flexural response is ap- 
parent long before the plate displacements reach the plate thick- 
ness indicating no significant membrane response involvement. 

The first plate/impactor separation occurs when the plate dis- 
placement is approximately equal to the plate thickness. Thus, 
even though membrane stretching is present, it is not the domi- 
nant response. 
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STRESS AND FAILURE PREDICTIONS 


Having characterized the displacement response of the 
clamped plate, a second analysis was performed in order to pre- 
dict stresses and damage within the plate due to the impact 
event. A new finite element mesh was developed since the size 
of the individual element in the mesh of figure 2 precluded re- 
alistic stress prediction. The finite element model utilized 
for stress and failure predictions is depicted in figure 6. The 
shaded elements in figures 6 were monitored for stress and damage. 

The refined model contains 521 unconstrained nodal points 
with three bending degrees of freedom per node. In addition, in- 
plane membrane degrees of freedom were retained in order to model 
any bending/extensional coupling introduced due to unsymmetric 
damage. Thus, the refined model contained 2605 active degrees of 
freedom. 

The inclusion of significantly more elements and degrees of 
freedom in the mesh introduced significantly higher flexural fre- 
quencies into the possible nodes of the model. As was stated 
previously, the solution algorithm utilized is sensitive to all 
possible dynamic nodes in the model. Thus, the maximum integra- 
tion time step was reduced in comparison to the model utilized 
for the displacement only model. In the analysis performed with 
the refined model, a time step no larger than 0.1 ysec was re- 
quired for numerical stability. 

Utilizing the finite element mesh of figure 6 and identical 
laminate and impact conditions as used previously, displacement, 
stress and damage were predicted for 100 ysec into the impact 
event. Using time step of 0.1 ysec, the analysis was performed 
for 1000 time increments. In the analysis, stresses were computed 
every 5 ysec or every 50 time steps. 

The lateral displacements of the impactor and node on the 
plate, which the contact spring is attached to, are plotted in 
figure 7. The impact mass is seen to be displacing along a 
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smooth path while the path of the plate is slightly irregular. 

This is similar to the type of response predicted with the coarse 
model. Comparing figure 7 and the first 100 ysec plotted in fig- 
ure 3, the two finite element meshes are seen to yield similar 
displacement-time responses. 

The displacement response of the plate through 100 ysec is 
depicted in figures 8 and 9. The displacements are plotted for 
lines through the plate center with figure 8 corresponding to the 
larger chord of the plate and figure 9 the shorter chord. Com- 
parison of figures 5 and 8 demonstrates that the coarse and re- 
fined finite element meshes yield the same maximum displacements 
for the plate. The curvatures produced in the refined mesh are 
significantly larger, however. 

The presence of high flexural nodes is easily discernable 
in figures 8 and 9. The plate response is seen to begin with a 
highly localized deformation at the point of impact. As the de- 
formation continues, it takes on the shape of a flexural wave 
spreading from the point of contact to the plate boundaries. As 
the wave passes through a point on the plate, the local curvatures 
reverse in sign. This is graphically demonstrated in figure 10 
where fiber direction stresses in the bottom 45° ply are depicted 
as a function of time. The stresses are plotted for element 79. 
The position in the nodel of element 79 can be found in figure 11 
where the finite element mesh is reproduced with element niimbers 
included. 

The stresses in figure 10 are seen to be compressive from 
the start of the impact event through approximately 33 ysec. The 
negative stresses during this period are due to the reversed cur- 
vatures which can be seen in figures 8 and 9. Beyond 33 ysec, 
the depicted stresses become positive and remain positive for the 
remainder of the 100 ysec analysis. At the point in time where 
the fiber direction stresses in the back ply of element 79 become 
zero, the same ply in adjacent elements have either positive or 
negative stresses depending upon their position relative to the 
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spreading flexural wave. At 33 ysec, the advancing flexural wave 
is located such that in element 79 the back 45° ply is unstressed 
in the direction of the fiber. Beyond 33 psec, the stresses in 
the back ply of element 79 become positive and oscillate somewhat 
due to the presence of high flexural modes and possibly induced 
damage. 

Damage which was predicted in the 100 ysec analysis is de- 
picted in figure 12. Only those elements which were monitored 
for stress and damage are included in figure 12. Thus, the ele- 
ments represented in figure 12 correspond to the shaded elements 
depicted in figure 6. Damage is first predicted at 35 ysec. At 
this time, the four elements surrounding the point of impact have 
experienced some damage. Elements 90 and 103 have experienced 
fiber failure in ply 1 (back face 45° ply) and matrix failure in 
ply 2 (0° ply) . Elements 91 and 102 have experienced matrix 
failure in ply 1 at 35 ysec. 

The occurrence of fiber failure in the back 45° ply of ele- 
ments 90 and 103 is interesting when compared with the fiber di- 
rection stresses of the back 45° ply of element 79 depicted in 
figure 10. Figure 10 indicates that at 35 ysec the ply 1 fiber 
direction stress in element 79 is less than 10 MPa while the same 
ply in elements 90 and 103 has failed in the fiber direction. The 
large difference in the stress state between these locations is 
due to the location of the spreading flexural wave. Since at 35 
ysec the transition from compression to tension has only just 
passed element 79, low stresses are to be expected. Elements 90 
and 103 are adjacent to the impact site where curvatures are max- 
imum and, hence, the stresses there can be expected to be large. 

Returning to figure 12 and the progression of impact induced 
damage, it is seen that at 40 ysec, elements 91 and 103 have ex- 
perienced matrix failure in ply 3. Ply 3 is oriented at -45° and 
is below the midsurface. Since ply 3 is oriented at 90° with re- 
spect to ply 1, the induced matrix failure suggests a similar 
loading state which promotes fiber failure in ply 1 at 35 ysec. 
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The combination of curvatures which produced large fiber direction 
stresses in ply 1 will produce stresses perpendicular to the fiber 
in ply 3 and, hence, matrix failure occurs in ply 3. 

The combined effects of compressive strengths being larger 
than tensile strengths and shifting of the neutral surface upward 
due to damage below the midplane, delays the onset of damage above 
the midsurface until 45 psec. At this time, elements 90 and 103 
experience fiber failure in ply 8 (top 45° ply) . The laminate re- 
mains in this damage state until 70 psec have elapsed in the impact 
event. 

At 70 psec into the impact event, matrix damage in ply 2 
spreads to elements 91 and 102. Ply 2 experiences additional ma- 
trix mode damage in elements 78 and 115 at 80 psec. At 90 psec, 
matrix mode damage occurs in ply 1 of elements 79 and 114. Final- 
ly, at 100 psec, ply 2 experiences additional matrix mode damage 
in elements 79 and 114. The progression of damage from the four 
elements surrounding the impact site (90, 91, 102, 103) to elements 
further away (78, 79, 114, 115) reflects the spreading of the flex- 
ural wave and the effects of included damage on subsequent plate 
response. The progression of damage in the 100 psec analysis is 
summarized in tables 3 and 4. Although interlaminar shear stresses 
were found to be as high as approximately 28 MPa, no delaminations 
were predicted during the 100 psec analysis. 


14 


DISCUSSION 


The impact analyses performed have generated information 
leading to a better understanding of the impact response of thin 
composite plates. The most interesting feature brought out in 
these analyses is the introduction of damage very early in the 
impact event. 

The displacement fields which were depicted in figures 8 
and 9 demonstrated that at times less than 100 ysec, the plate 
responded with very high local curvatures at the point of impact. 
These large curvatures produced considerable local damage. As 
the impact event progressed, the displaced shapes of the plate 
began to approach what appeared to be the fundamental flexural 
mode of the clamped-clamped plate. The damage had appeared be- 
fore the fundamental mode could be attained, however. Thus, when 
the plate response degenerates into this mode, the properties of 
the plate have changed from its virgin state. This indicates that 
the first mode response of the plate must be different in the im- 
pact event than it would be under static or quasi-static loading 
since static static loads cannot produce the very high local 
curvatures and resultant damage. 

The prediction of damage in the current analyses is con- 
siderably different than that experienced in the first phase of 
this effort (ref. 14). In both models, the planar area of dam- 
age was approximately 1.0 cm x 1.5 cm. In the earlier effort, 
the analyses predicted such extensive damage at the point of im- 
pact that the elemental stiffness matrices became singular. Ad- 
ditionally, the damage occurred earlier such that the maximum 
time in the solution was 75 ysec. The current analysis produced 
much less damage through the plate thickness. 

The damage predicted in the current effort more closely 
simulates the damage produced experimentally in tests performed 
at NASA Langley than did the earlier effort. The impact parame- 
ters and laminate configuration used in the analyses were pre- 
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scribed to analyze experimental work from which qualitative dam- 
age descriptions were available. The experimental work produced 
matrix damage and delaminations but not significant fiber damage. 
Thus, the earlier model severely over predicted the extent of 
damage. The current analysis predicted much less extensive dam- 
age and, hence, more closely predicted the experimental work. 

The primary reason for the reduction in the extent of pre- 
dicted damage relates to the way in which the impact load is in- 
troduced into the plate. In the earlier effort, the impact mass 
was lumped at nodes of the plate and given an initial velocity. 
Thus, a step loading was applied. In the current model, a non- 
linear contact spring is utilized. This allows for a more real- 
istic gradual introduction of load. The contact spring stiffness 
is initially zero but increases with increasing displacement of 
the impactor relative to the plate. 

Another factor which must be considered is the relative sizes 
of the individual elements surrounding the impact site. Although 
the two models contained similar numbers of flexural degrees of 
freedom (1953 earlier effort, 1563 current effort) , the current 
element requires midsize nodes. Therefore, fewer elements can be 
included for the same nximber of degrees of freedom. Thus, in the 
earlier effort, the elements at the impact site were approximately 
0.16 cm X 0.16 cm while in the current effort the impact site el- 
ements were approximately 0.42 cm x 0.42 cm. Therefore, some of 
the local stress response may have been lost in the current ef- 
fort. It is believed that element size caused the lack of inter- 
laminar damage in the current effort. 

The obvious solution to this problem would be to increase 
the nvimber of finite elements in the model and decrease the size 
of the elements local to the impact site. The difficulty is the 
magnitude of the finite element model that this would induce. 

In the current model, there were 2605 active degrees of freedom 
(including in-plane displacements) yielding 5210 simultaneous 


16 


differential equations (see Appendix C) . This model approached 
the maximum size which was pratical to work with. 

The removal of in-plane degrees of freedom would signifi- 
cantly reduce the size of the problem while eliminating any bend- 
ing/extensional coupling effects. This could be accomplished in 
a clamped plate by reformulating the laminate bending properties 
to enforce bending about the neutral rather than midsurface. Such 
a modification would be more difficult to define for a simply sup- 
ported plate. 

The size of the current solution was one reason for terminat- 
ing the analysis at 100 ysec. As was stated, this corresponded to 
1000 time increments. The current analysis is very costly to ex- 
ecute which must be considered typical of transient dynamics anal- 
yses. In addition to the time and expense of operating the code, 
the displacements at 100 ysec were approaching 75% of the plate 
thickness. The validity of the solution becomes suspect as the 
displacements reach the plate thickness and membrane stretching 
effects come into play. For the eight ply laminate utilized in 
the analysis, it can be expected that membrane effects will be- 
come significant if not dominant at times later than 100 ysec. 
Since the analysis cannot handle large deflection effects, it was 
deemed prudent to stop the analysis at 100 ysec. 
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RECOMMENDATIONS 


The impact analysis methodology which has been incorporated 
into the CLIP II computer code has identified some significant 
phenomena regarding low velocity impact of thin composite lami- 
nates. Specifically, significant damage producing curvatures 
have been identified at times very early in the impact event. 

This type of deformation, which is characteristic of very high 
dynamic flexural mode exitation, has been shown to promote damage 
within the first 100 ysec of a typical impact event. The trans- 
ient dynamics finite element analysis has proven to be costly, 
however, and therefore future efforts should be directed towards 
defining realistic simplifications to the analysis in order to 
obtain a practical design/analysis tool. 

The results obtained utilizing the CLIP II code can be seen 
to indicate two distinct modes of plate response and, hence, dam- 
age producing mechanisms. In the very early stages of an impact 
event, very high frequency flexural effects dominate the response. 
At considerably later times, the plate response degenerates into 
what appears to be a fundamental flexural mode response, possibly 
coupled with large deflection effects. The significance of the 
large deflection membrane response will be determined both by the 
characteristics of the plate (i.e. thickness) and the impactor en- 
ergy and velocity. Simplifications to the impact analysis will 
necessarily be directed towards these two regimes independently. 

The early response of a thin plate subjected impact loading 
is characterized by high local curvatures. This type of response 
should be present in impact of thicker plates although perhaps not 
as pronounced. Simplifications of early time response analysis 
should take advantage of the localized nature of the high curva- 
tures. Early in the impact event, significant portions of the 
plate are unaffected. As the event progresses, the portion of the 
plate which is unaffected decreases until finally the entire plate 
is deforming. It should be possible to define characteristic di- 
mensions of the plate as functions of time corresponding to the 
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deformed regions of the plate. This information could then be 
used to reduce the size of the required finite element model. 
Alternately, this data could be utilized to develop one-dimensional 
models where a first mode solution would be used to model the de- 
formations of the affected regions of the plate. Either of these 
approaches should greatly reduce the time and cost of performing 
the analyses for early time response. 

The later time response might be simulated using a first mode 
flexural analysis with large deflection capabilities. The models 
would need to include the effects of any damage which occurred 
during the early times of the impact event. This approach would 
require that late time plate response does, in fact, degenerate 
into the fundamental mode of the plate structure. 

Simplifications of either short or long term impact response 
will require considerable study of impact phenomena. Realms of 
applicability of possible simplifications must be defined in or- 
der that the analyses be realistic. The information to be gath- 
ered could be obtained experimentally or analytically. Experimen- 
tal efforts would need to monitor the plate response in such a way 
as to yield both displacement vs. time and displacement vs. posi- 
tion. Analytical efforts would also have to characterize the time- 
position characteristics of the displacement response. 

The types of analysis simplifications described would allow 
for the development of cost effective design/analysis tools which 
could be utilized to investigate impact events. The current me- 
thodology must be considered too costly to be utilized routinely. 
Simplified analytical methodologies with known limitations would 
be a valuable contribution in the continuing study of low velocity 
impact. 
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Table 1. Ply Properties (T300/5208, = . 
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2 

12 

12 

13 

23 


153.0 GPa 
10.9 GPa 

.3 

5.6 GPa 
5.6 GPa 
4.2 GPa 
1.55 g/cc 
.0132 cm 


Ply Strengths 


^IC 
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MPa 

®1T 
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MPa 
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= 

96 

.5 MPa 

^2T 

= 

27 

.6 MPa 

^12 

= 

62 

MPa 

^13 

= 

62 

MPa 

^23 

= 

62 

MPa 

®IF 

= 

62 

MPa 
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Table 2. Impactor and Contact Spring Parameters 


Impactor Mass = 16.45 g (1.59 cm dia. steel ball) 

Initial Velocity = -9.4 m/s 
Initial Position = 0.0 


Spring constant (K) = 6.55 x 10^ (N/m^) 

Loading exponent (n) = 1.5 
Unloading exponent (m) = 2.5 

4 1-i 

Indentation constant (c) = 6.067 x 10 (m ) 
Indentation exponent (i) =2.4 
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Figure 5. Lateral Displacement of the Major Chord 
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Figure 8. Lateral Deflection of the Major Chord 
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APPENDIX A 


THE ISOPARAMETRIC THICK SHELL ELEMENT 


The element used in CLIP II is an extension of the eight noded 
isoparametric plate bending element with shear deformation devel- 
oped by Hinton et. al (ref. 19) . Membrane effects and bending- 
extensional coupling have been added to the basic element presented 
in reference 19, resulting in an element well suited for the study 
of composite plates. 

Although the Kirchoff hypothesis that lines originally normal 
to the plate remain normal after deformation is not used, other as- 
sumptions have been made. They are: 

1. Lateral displacements are small. 

2. Lines originally normal to the plate remain straight. 

3. Normal strains and stresses are negligible. 

Shear deformation has been included in this element by allow- 
ing the rotations of the normal about the global x and y axes, and 
the lateral plate displacement to very independently. By assuming 
a linear variation of u and v with respect to the z axis (assump- 
tion 2) , the element may be integrated explicitly through the thick- 
ness, reducing the three dimensional problem to two dimensions. 

While it is known that normals do not stay straight during de- 
formation, the rotations about the coordinate axes may be thought 
of as average values, and the actual nonuniform through- the-thick- 
ness shear distribution accounted for by an alternate procedure. 

The displacement field for this element is given by: 


u(x,y,z) = u (x,y)+z0 (x,y) 

O A 

(A.l) 

v(x,y,z) = v^(x,y) -z9y (x,y) 

(A. 2) 

w(x,y,z) = w^(x,y) 

(A. 3) 
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where u^, and are the displacements of the mid-plane, and 
0^ and 0y are the rotations in the x-z and y-z planes, as shown 
in figures A.l and A. 2. 

Using this displacement field and assumptions 1 and 3, the 
strains are given by 



These strains may be separated into two parts, one accounting 
for membrane strains and the other for bending strains. 

{T} = + {"Eg} = [Bj^] {u} + [Bg] {u} (A. 5) 

where 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

0 0 0 
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0 0 0 


0 +Z ^ 
3y 


0 ^ 


ay 


0 

0 -zf 

3y 

3x 
0 

-1 


+1 


-z 


(A. 7) 


and 


{u}'^ = {u^(x,y) v^(x,y) w^(x,y) 0^(x,y) 6y(x,y)} (A. 8) 


The displacement field is approximated with parabolic shape 
functions, using the displacements and rotations at the eight nodes 
shown in figure A. 2. To aid in itegrating over the area of the 
element, the original element geometry is mapped onto a two unit 
square in an r-s coordinate system, as shown in figure A. 2. Be- 
ing isoparametric, this coordinate mapping is done with the same 
functions used to approximate the displacement field. These func- 
tions are of the serindipity type, and have the form 


h. 

1 


• (1+r^r) • (1+Sj^s) • (r^r+Sj^s-1) 


i=l,2,3,4 


(A. 9) 


and 


hi = -| • (l+rfr) • (l+SiS)-(l-ri^r^) • (l-Sf^s^) i=5,6,7,8 (A. 10) 

Using these shape functions, any function of r and s may be 
approximated by 

8 

f(r,s) = Z h.f. (A. 11) 

i=l ^ ^ 
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Further, partial derivatives with respect to r and s are given by 


M 

3r 


8 

Z 

i=l 


3r 


f . 

1 


3s 


8 

E 

i=l 


3h. 
1 

3s 


f . 
1 


(A. 12) 


(A. 13) 


Since strains have been defined in terms of the derivatives 
with respect to x and y, a coordinate transformation for the deriv- 
atives is needed. This transformation is found by applying the 
chain rule for derivatives. 



" 3x 

3y‘ 

/ 3 I 

1 

3 

3r 

3r 

1 ^ 1 

3x 



r 


3x 


1 i_ i 

3 

_ 

3s _ 

1 / 

3y 


(A. 14) 


The Jacobian transformation matrix [J] is evaluated with the 
aid of (A. 12) and (A. 13). 


[ J] = 


■ 3h^ 

31.2 

3h 

3r 

3r . 


9h^ 

9h2 

9h 

L 3s 

3s 

Ts 


xi y^ 


X2 Y2 


• • 


^8 ^8 


(A. 15) 
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Partial derivatives with respect to x and y are then given 
by inverting (A. 14) and substituting (A. 12) and (A. 13). 


/ ^ ' 


■ 3hj^ 


^^81 


f \ 

f 

1 8x 

-1 

3r 

_ • • • « 

3r 

3r 


^1 

1 9f 

■ = [J] 

3h^ 

3h2 

00 


f 

(by 

1 

L 3s 

JT~ 

3s _ 

< 

2 


= [p]|f 1 


(A. 16) 


V ^8/ 

The strain-displacement relationship for the element is now 
obtained by substituting (A. 16) into (A. 6) and (A. 7). 


Iil= ([B„J + [Bgl) jul 


(A. 17) 


with 


l“( = Is S % 


0.. e„ 0 


°8 °8 *^8 ^8 ^8 ) (A. 18) 


! (A.: 


= 


Pii ° 


0 0 0 Pj^2 ° 


P 2 I 


^21 ^11 


0 0 0 0 


■22 


0 0 0 P 22 P 12 


0 0 0 0 


0 0 0 
0 0 0 
0 0 0 
0 0 0 


18 


0 0 0 0 

0 Pjg 000 

P 28 P18 0 ° ° 

0 0 0 0 0 


(A. 19) 
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and 




0 0 0 
0 0 0 0 
0 0 0 zp. 

0 0 Pii +1 
.0 0 0 


zpj^j^ 0 0 0 0 zp 


12 


0 0 0 0 


-zp 


-ZP21 


-zp, 0 0 0 zp„ -zp 


22 


28 


21 ‘‘'11 


22 ‘*'12 • 


0 0 Pj^2 +1 


-1 0 0 


0 -1 


t-22 


..000 zpjg 0 

.. 0 0 0 0 -zp 

•000 ZP28 -ZP18 

. . 0 0 Pj^8 +1 0 

. 0 0 Pj8 0 -1 J 


(A. 20) 


where is the element of the matrix [P] defined by (A. 16) 

Using standard finite element techniques, see references 20 
and 21, the stiffness matrix for this element may be calculated 
from 


[K] = /^°^[B]'^[D] [B]dv 


(A. 21) 


where [D] is the 5x5 stress-strain relationship defined by 


{a} = [D]{e} 


(A. 22) 


where 


ip 

{a} = {o a T T T } 

X y xy xz yz 


(A. 23) 


and 


{e} ^ ^x^y''^xy"*^xz"'^yz^ 


(A. 24) 


Since the matrix [B] is the sum of two matrices, the integral 
in (A. 21) may be rewritten as 


[K] = /"'°^[B^]'^[D] [Bj^ldv + + 

[Bj^ldv + /’^°^[Bg]'^[D] [Bgldv 


(A. 25) 
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The total element stiffness matrix may then be thought of as the 
sum of four stiffness matrices 

[K] = [K^] + [K 2 ] + [K 3 ] + [K^] (A. 26) 

considering the first of these four matrices 

[K^] = [B^ldv (A. 27) 

From (A. 19), it can be seen that the matrix [B^^] is independent of 
z, allowing (A. 26) to be rewritten as 

[K^] = [Bj^]dA (A. 28) 

This matrix is recognized as the plane stress stiffness matrix, 
and includes only in-plane effects. The upper 3x3 of [D] is the 
only portion of this matrix that is utilized, and is recognized 
as the [A] matrix from laminated plate theory, when multiplied by 

t. 

The metrics [K 2 ] and [K^] are considered together, since it 
can be shown that 

[K 2 ] = [K 3 ]'^ (A. 29) 

The expression for [B„] , given by (A. 20), may be expressed as the 

13 

product of two matrices 

[Bg] = [Z][B*] (A. 30) 

where 
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z 0 0 0 0 

0 z 0 0 0 


[Z] = 


0 


0 z 0 


0 0 0 1 
0 0 0 0 
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0 
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0 
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0 

?22 

0 
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° P18 
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... 0 

0 

0 0 

"P28 

. . . 0 

0 

° P28 

“P18 

. . . 0 

0 

P18 ^ 

0 

. . . 0 

0 

P28 ° 

-1 _ 


(A. 31) 


(A. 32) 


This allows the expression for [K 2 ] to be written as 
[K^] = [D]dz) [Bj^ldA 


(A. 33) 


The innermost integral of (A. 33) may be evaluated explicitly, and 
represents bending-extensional coupling. The upper 3x3 of the re- 
sulting integration is the [B] matrix from laminated plate theory. 
For homogeneous isotropic materials or balanced symmetric laminates, 
this integral will vanish. 

The expression for [B„] given by (A. 30) is also utilized to 

n 

evaluate [K^] . Substituting this expresison yields 

[K^] = [Z]dz) [B*]dA (A. 34) 
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Again, the innermost integral of (A. 34) may be evaluated expli- 
citly. This integral represents the bending properties of the 
element, and may give rise to bending- twisting coupling, depending 
on the nature of [D] . The upper 3x3 of this matrix is the [D] ma- 
trix from laminated plate theory, the lower right 2x2 submatrix 
of [D] contains the shear terms. It is this submatrix which must 
be corrected to account for nonuniform shear distributions in the 
z direction. For homogeneous isotropic plates, the shear distri- 
bution is known to be parabolic, resulting in a correction factor 
of 5Gt/6. For laminates, however, the shear distribution is not 
so well defined, requiring that more complicated procedures be used. 

The procedure used in CLIP II is based on the method developed 
by Cohen in reference 19. Cohen presents a method for determining 
the constitutive equation for transverse shear without specifying 
a priori the through-the-thickness distribution. This method in- 
volves the application of Castigliano' s theorem to minimize the 
shear strain energy, which is defined in terms of a Taylor series 
containing arbitrary constants. The details of this procedure are 
beyond the scope of this report, but can be found in reference 19. 
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Figure 
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Figure A. 2. Extensional Strain in the Y Direction 
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APPENDIX B 

THE NONLINEAR CONTACT SPRING 


CLIP II models the impactor as a lumped mass which is tied 
to the plate at a single node by a nonlinear spring. This spring 
contains hysterisis, and has the ability to vanish if the plate 
and impactor have separated. The general force-displacement re- 
lationship for this spring were determined by Sun and Yang (refs. 
16 and 17) and are 

F = K6^ (B.l) 


if the spring is loading, and by 


F 


F 

max 


6-6 


m 


'6 -6 ' 

max o 


(B.2) 


if the spring is unloading. The quantity 6^ is the permanent de- 
formation and is determined as a two parameter curve fit from ex- 
perimental data. The form used is 


6 = C6^ 

o max 


(B.3) 


This quantity represents how deeply the impactor has indented the 
plate. The quantities K, n, m, C and i must be found experimentally. 
Yang and Sun, (refs. 16 and 17) , have shown that this formulation 
gives reasonable results for composite plates being impacted by a 
spherical impactor. If 6 becomes less than 6^, the plate and im- 
oactor have separated, and the contact force is set to zero. 

A typical force-deflection diagram is shown in figure B.l. 

This curve shows a spring that is loading from point A to B. At 
point B, the plate has undergone a permanent indentation 
so that it will unload along curve 2 to point C. At point C, the 
plate and impactor have separated, allowing further separation to 
follow curve 3 to point D. Curve 4 shows the gap between the plate 
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and impactor closing until contact is again made at point C. • The 
spring then loads along curve 5 until it reaches point B. At B, 
it continues loading along the original curve until point E, where 
the spring begins to unload again. The plate now has a new perma- 
nent indentation causing the spring to unload along curve 7 

to point F, where it reverses and begins to load again. During 
the course of an impact event, this loading, unloading and separa- 
tion may occur many times. 
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APPENDIX C 

NUMERICAL SOLUTION OF THE TRANSIENT DYNAMIC PROBLEM 

The general form of the equations of motion for a discretized 
elastic system may be written as 

[M]{x} + [C]{x} + [K]{x} = {F} (C.l) 


where 


[M] 

= 

mass matrix 

[C] 

= 

damping matrix 

[K] 

= 

stiffness matrix 

{x} 

= 

acceleration vector ( 

{k} 

= 

velocity vector (^) 

{x} 

= 

displacement vector 

{F} 

= 

applied load vector 


This equation represents a set of n simultaneous second order 
differential equations. This set of n second order equations is 
most easily solved numerically, if it is rewritten as a set of 2n 
first order equations. This may be done by defining 

{x} = {y} (C.2) 

Substituting (C.2) into (C.l) and solving for {y} yields 

{y} = [M]'^({F}-[C] {y} - [K]{x}) (C.3) 

Two assumptions have been made to ease the computational ef- 
fort. First, the mass matrix is assumed to be diagonal (Itimped 
mass) , and secondly, Rayleigh damping is used where the damping 
matrix is assumed to be related to the mass and stiffness matrices 
by 
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[C] = a[M] + 3[K] 


(C.4) 


Where the constants a and 6 are related to the material and struc- 
tural damping of the system. Incorporating these assumptions into 
equation (C.3) and combining the result with equation (C.2) gives 


« 

X 


{<i>} 


[4)] [I] 


X 

y 




[M]“^[K] (a+3[M]“^[K] ) 


y 


(C.5) 


Since the global stiffness matrix is merely the sum of the 
individual element stiffness matrices, equation (C.5) may be eval- 
uated by summing the contribution due to each element without ever 
assembling the global stiffness matrix. The global mass matrix 
must, however, be assembled. Inversion of the mass matrix is sim- 
plified due to the fact that it is diagonal. Internally, CLIP II 
treats the mass matrix as a vector, which reduces both the storage 
requirements, and the number of multiplications needed for a matrix 
multiplication. A computational simple expression for y^ may then 
be written as 


1 n 

"ii ^ j=i 


K..X. - aM..y. - 3 ? k..y.) 
ID 3 11 1 ID D 


(C.6) 


Equation (C.5) is now in a form which allows a suitable numer- 
ical integration scheme to be used to step out the dynamic solution 
in time. The scheme used in CLIP II is the fourth order Runge- 
Kutta-Gill method, and is implemented in the subfoutine RKGINT. 

This integration formula predicts the solution at time i+1 wing 


= {u}. + f [{K^} + 2(1^) {K^} + 2 (H^){K 3 ) 

(C.7) 

+ {K4}] 
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where 


{u}^ = {x y} 

{K^} = {g (t^, {u}^) } 

{K^} = {g(t^+|At, {u}^4^At{K^}} 

{K^} = {g(t^-f^At, {u}^+(--|4^)At{K^}+(l-^)At{K2}} 
{K^} = {g{t^+At, {u}^~At{K2}+(l-^^)At{K2}} 


and {g(t^, {u}^) } represents equation (C.5) evaluated at time t^ 
with the solution vector {u}^. This evaluation is performed within 
CLIP II by the subroutine DERIVS. This scheme requires that the 
initial state of the system, both velocities and displacements, be 
specified. 

The time increment. At, may vary from step to step, as only 
the last solution is required to predict the next solution. If, 
however, the time increment is relatively large, it may be neces- 
sary to periodically correct the predicted solution. This is done 
in CLIP II by the subroutine CORECT, which incorporates Hamming's 
fourth order correcting algorithm. This algorithm is of the form 


{u} 


i+1 




(C.8) 


Notice that equation (C.8) is an implicit formula which re- 
quires an iterative solution. Further, it is dependent upon the 
last three solutions, which requires the time increment to be con- 
stant for at least three time steps before the correcting is per- 
formed. This procedure is very time cons\iming, and should only be 
done if absolutely necessary. Reference 20 provides a good dis- 
cussion on both the R-K-G predictor and Hamming's corrector. 
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APPENDIX D 

LAMINATE FAILURE CRITERIA 


The failure criteria utilized in CLIP II to predict damage 
are based upon Kashin's failure criteria for unidirectional fiber 
composites (ref. 23) . The criteria allows for the prediction of 
both the presence of failure and identification of the mode of 
failure. Modes of failure which can be predicted include fiber 
breakage, ply splitting parallel to the fibers, and delaminations 
between plies. 

The stresses predicted within CLIP II are depicted in figure 
D.l. The 1-direction coressponds to the fiber direction in a ply 
and ^23 ®'^^®sses are transverse shear stresses. The fail- 

ure criteria make use of eight material strength parameters which 
are: 


SfT 

= 

fiber 

direction tensile strength 


^IC 


fiber 

direction compressive strength 


^2T 


tensile strength perpendicular to the 

fibers 

®2C 

= 

compressive strength perpendicular to 

the fibers 

^12 

s 

shear 

strength in the 1-2 plane 


^13 

= 

shear 

strength in the 1-3 plane 


^23 

= 

shear 

strength in the 2-3 plane 


®IF 

= 

interlaminar shear strength 



The failure criteria are separated into the three different 
modes. For fiber breakage, the criteria are 



(D.l) 


at failure for > 0, and 
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at failure for a^^<0. 

For predicting splitting parallel to the fibers, the criteria 

are 



at failure for 022 > 0/ 



(D.3) 


(D.4) 


at failure for 022 < 0. 

Delaminations are predicted utilizing a simple quadratic in- 
teraction formula: 


o 


2 

13 


+ o 


23 


S 


2 

IF 


1 


(D.5) 


at failure. 

In the prediction of failure, equations (D.l) - (D.4) are 
evaluated at the top, middle and bottom of each ply in the laminate 
and equation (D.5) is evaluated at both the top and bottom of each 
ply. 
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APPENDIX E 

INCORPORATION OF DELAMINATION EFFECTS 


The methodology adopted for incorporating the effects of dam- 
age in the impact analysis involves modification of the elemental 
stiffness matrices. For the case of ply damage, the reduction in 
stiffness is accounted for simply by specifying appropriate ply 
moduli as zero and reformulating the laminate constitutive matrices. 

Elemental modifications to account for delaminations are not 
so easily defined since the result of a delamination is a change 
in both the bending and transverse shear stress distributions. The 
change in bending stress is a function of the presence of shear and 
hence the two are coupled. A delamination does not affect the 
stiffness of a bending element in the presence of pure bending. 

Consider a homogeneous beam with a central delamination sub- 
jected to pure bending (no shear) . Since there are no shear 
stresses, the plane of the delamination is stress free. Thus, 
the delamination can have no effect on the stress field or bending 
rigidity. 

If a homogeneous beam is subjected to combined bending and 
shear, the resulting stresses are those depicted in figure E.2. 

The parabolic shear stress distribution is characteristic of the 
linear bending stress distribution. When a delamination is intro- 
duced, the bending and shear stresses must change, primarily be- 
cause the shear stresses must vanish at the delamination. 

The bending stresses in the delaminated beam subjected to both 
bending and shear are shown in figure E.3. In addition to the ex- 
pected bending stress variation, singular forces at the ends of 
the delamination are present which balance the difference in total 
force at each end of the delaminated beam. Thus, horizontal equi- 
libri\im is maintained. Since the magnitude of the singular forces 
balances the force imbalance in bending, the singular force is 
simply the shear force which existed on the plane of the delamina- 
tion prior to the delamination. 
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The stresses depicted in figure E.3 can be separated into two 
parts consisting of a normal bending stress distribution and a mem- 
brane stress distribution. This has been done in figure E.4. The 
magnitude of the singular forces is now seen to be simply the dif- 
ference in membrane force along the length of the partial beams. 

The moments depicted in figure E.4 must be balanced by the 
transverse shear forces. Without consideration of the singular 
forces, the moment gradients depicted are not in balance with the 
shear forces. This is remedied by distributing the singular forces 
over the partial beams as depicted in figure E.5. The singular 
forces have been moved to the midsurfaces of the two partial beams 
and moments added to account for the change in the line of action 
of the moved forces. The result of the moments added serves to 
increase the moment gradients across the length of the partial 
beams. This occurs since the smaller end moments have been reduced 
while the larger end moments have been increased. It can be shown 
that the moment gradients are now equal to the shear force and 
hence the model is in equilibrium. The resulting bending stress 
distributions are depicted in figure E.6. Shear stresses are also 
shown in figure E.6. The peak shear stresses are identical to 
those shown in figure E.2. 

The model described for the effect of bending and shear in 
the presence of a delamination has been utilized to define proper- 
ties of a composite laminate with a delamination. The development 
of the analysis for a laminate requires the partial laminate con- 
stitutive relations and expressions for transverse shear stresses 
as a function of the average rotation. 

Utilizing the work of Cohen (ref. 19) , the relation between 
interlaminar shear stress and transverse shear rotation is de- 
fined as 


T 

XZ 

^ i 

Y 

'xz 

T 

zy 

II 

1 — 1 
w 
1 1 

1 

II 

Y 

yz 
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' X z 

where } are average shear rotations and the superscript "i" 

'yz 

identifies the location through the thickness. Using (E.l), sin- 
gular forces due to a delamination at interface "i" can be deter- 
mined as 


F ® = /^ T dA 
X A xz 


and 


F = f, T dA 
y A xz 


or 


S V 


X 


[R]"" /, 


/ 1 


xz I 


Y. 


dA 


yz 


(E.2) 


In the CLIP II code, stresses are computed as an average over the 
element so that equations (E.2) reduce to 




pxz 

I ^yz 


(E.3) 


where A® is the planar area of the element. 

The forces defined in (E.3) are converted into stress and mo- 
ment resultants in a fashion similar to that shown in figure E.4. 
Stress resultants become 


N ° 

X 


1/Ay Fg 

N ^ 

XT 


1/Ax F 
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where Ax and Ay are the elemental widths in the x and y directions. 
The requirement that elements to be monitored for stress be defined 
as shown in figure 1 stems from the use of equations (E.4) and the 
need to define Ax and Ay for an element. 

Moment resultants are defined as 


M 

X 

= (z^-z.) 

N 

X 

M ® 

\T 

S 1 

N 

XT 


where Z. is the location of the delamination and Z is the mid- 
1 s 

surface of the partial laminate either above or below the delamina- 
tion (fig. E. 7) . 

The stress resultants (E.4) and moment resultants (E.5) are 
applied, half at each end, to the partial laminates with a linear 
variation, plus to minus, similarly to that shown in figure E.4. 
Applying half at each end implies that the singular forces are dis- 
tributed equally at each end of the delamination. 

Using the stress and moment resultants, strains and curvatures 
are defined using constitutive relations for the partial laminates. 
Thus 
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A 

B 

-1 

N ^ 

X 

— 




X 

s 

K 


B 

D 


M ® 

X 





X 


These strains and curvatures are tllen utilized in defining the 
variation of strain through the thickness of the total laminate. 
The strain distribution through a plate with a delamination takes 
the form 


e = e + Zk + + (Z-Z )k^ 

o s 


(E.7) 
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where and may be different for the various sublaminates-. 

Thus, the strains are piecewise linear through the thickness of a 
delaminated plate element. 

This poses a fundamental problem in the plate element since 
different curvatures are allowed at different locations through 
the thickness and because jump discontinuities are allowed in 
strain through the thickness. Neither of these phenomena can be 
handled with a single plate element modelling all of the sublami- 
nates. To model these effects properly requires a separate plate 
element for each sublaminate. This approach leads to modelling 
each ply in the laminate as a separate plate element which then 
becomes such a large finite element model as to be deemed impossi- 
ble to utilize. 

The approach taken has been to include the strains and cur- 
vatures due to the singular forces while maintaining a single value 
for total rotation and midplane deformation. This leads to an ef- 
fective reduction in the total rigidity of an element thus incor- 
porating the effects of a delamination. 
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Figure E.3. Bending and Shear with a Delamination 


66 



( 

B t 

B : 

Figure E. 








I 

I 


. Decomposition of Bending Stresses with 
Shear and a Delamination 
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Figure E.5. Distributing Singular Shear Forces 




Figure E.6. Bending and Shear Stresses with a 
Delamination 
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APPENDIX F 


CLIP II PROGRAM USER'S GUIDE 


PROGRAM CONTROL DATA 


Card 1 


Title Card 


20A4 


Columns 



Contents 

1-80 

TITLE (20) 


80 character alphanvimeric title 

Card 2 


Save & 

Restart Flags 215A 

Coliomns 



Contents 

1-5 

IREST 


If IREST=0, this is a new analysis 




If IREST=1, restart the analysis 
from FILE08 

6-10 

I SAVE 


If ISAVE=0 , do not generate restart 
file (FILE08) 




If ISAVE=1, generate a restart file 

Card 3 


Step 

Control Data 315 

Columns 



Contents 

1-5 

NTSTEP 


The number of time steps to be taken 

6-10 

NPRINT 

■ 

The number of time steps pe^j^solu- 
tion printout (every NPRINT^ time 
step will be printed) 

11-15 

NSTRES 


The number of time steps per stress 
analysis (a stress analysis will be 
performed after NSTRES time steps) 

Card 4 


Printout Control Flags 415 

Columns 



Contents 

1-5 

LAMPRT 


Lamina stress-strain matrix print- 
out flag 

6-10 

VELPRT 


Velocity printout flag 

11-15 

ACCPRT 


Acceleration printout flag 

16-20 

STRPRT 


Detailed stress printout flag 


For any of the above flags, a value 
of 1 turns on the printout option, 
and a value of 0 suppresses the 
the printout 
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I. PROGRAM CONTROL DATA (continued) 


Card 5 


Time Step Data 

2F10. 0 

Columns 


Contents 


1-10 

TIME 

The initial time 


11-20 

DELT 

The time increment 


Card 6 

Correcting Algorithm Control Data 

215, FIO.O 

Col\imns 


Contents 


1-5 

NCOREC 

The number of time 

steps per cor- 



recting step. If NCOREC=0 , no 



correcting will be 

performed 

6-10 

NITMAX 

The maximum nvimber 

of iterations 



per correcting step 

11-20 

CONV 

The corrector convergence criterion 

Card 7 


Stiffness Integration Control 


Columns 


Contents 



1-5 lORDER The order of the nvimerical integra- 

tion used to evaluate the element 
stiffness matrices, may be either 
2 or 3 


II. MATERIAL DESCRIPTION INPUT 


Card 1 


Number of Materials 


15 


Columns Contents 

1-5 NMAT The nvimber of different materials 

used in the laminate. 0<NMAT<5 


72 


II. MATERIAL DESCRIPTION INPUT (continued) 


Card 2 Material Properties for Material 1 7F10. 0 

Columns Contents 


1-10 

El 

Longitudinal Young's 

modulus 

11-20 

E2 

Transverse Young's modulus 

21-30 

NU12 

Poisson's ratio in the 1-2 plane 
(load-strain) 

31-40 

G12 

Shear modulus in the 

1-2 plane 

41-50 

G13 

Shear modulus in the 

1-3 plane 

51-60 

G23 

Shear modulus in the 

2-3 plane 

61-70 

RHO 

Density 


Card 3 


Strengths of Material 1 

8F10. 0 

Columns 


Contents 


1-10 

SIC 

Longitudinal compressive strength 

11-20 

SIT 

Longitudinal tensile 

strength 

21-30 

S2C 

Transverse compressive strength 

31-40 

S2T 

Transverse tensile strength 

41-50 

S12 

Shear strength in the 

1-2 plane 

51-60 

S13 

Shear strength in the 

1-3 plane 

61-70 

S23 

Shear strength in the 

2-3 plane 

71-80 

SIF 

Interfacial shear strength 


If more than one material is used, cards 2 and 3 are repeated 
NMAT times. 


III. LAMINATE DESCRIPTION INPUT 


Card 1 


Number of Plies 

11 

Colvimns 


Contents 


1-5 

NPLY 

The number of plies in 
inate. 0<NPLY<25 

the lam' 
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III. LAMINATE DESCRIPTION INPUT (continued) 


Card 2 


Ply Description 


2F10.0, 15 

Columns 


Contents 



1-10 

THICK (I) 

The thickness of ply 1 



11-20 

THETA (I) 

The orientation of ply 

1, 

measured 



counterclockwise from 

the 

positive 



X axis 



21-25 

MID 

The I.D. number of the 

material used 



for ply 1 



Card 

2 is repeated NPLY times. 




IV. DAMPING INPUT 

Damping Parameters 2F10. 0 

Contents 

The mass matrix multiplier 
The stiffness matrix multiplier 


Card 1 
Columns 

1-10 ALPHA 

11-20 BETA 


V. NODAL COORDINATE INPUT 


Card 1 


Number 

of Nodes 15 

Col\imns 



Contents 

1-5 

NN 


The number of nodes used in the model 

Card 2 


Coordinate Input 15, 2F10.0 

Columns 



Contents 

1-5 

N 


Node number. N is only used for clar- 
ity when looking at the coordinate 
input. The M^h coordinate card will 
be used for the coordinates of node 
M regardless of the value of N. 

6-15 

X 


The X coordinate of node N 

16-25 

y 


The Y coordinate of node N 

Card 

2 will be 

repeated NN 

times. 
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VI 


CONSTRAINED DEGREE OF FREEDOM INPUT 


Card 1 Niomber of Constrained Nodes 15 

Columns Contents 


1-5 NCDF 


The number of nodes constrained in 
direction 1 (global x direction) 


Card 2 
Columns 

1-5 ICDF(l) 

6-10 ICDF(2) 


Constrained Nodes 1515 

Contents 

1st node constrained in direction 1 
2nd node constrained in direction 2 


71-75 ICDF(15) 15th node constrained in direction 3 

Card 2 is repeated until NCDF nodes have been specified. 

Card sets 1 and 2 are repeated for directions 2 , 3 , 4 and 5. 

Direction Degree of Freedom 


1 

2 

3 

4 

5 



VII. ELEMENT DEFINITION INPUT 


Card 1 
Columns 
1-5 NE 


Number of Elements 15 

Contents 

The number of elements used in the 
model 
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VII. ELEMENT DEFINITION INPUT (continued) 


Card 2 


Element Definition 

915 

Columns 



Contents 


1-5 

NI 

1st node 

in the element 

definition 

6-10 

NJ 

2nd node 

in the element 

definition 


36-40 NP 8th node in the element definition 

41-45 MON Stress analysis monitor. If MON=0, 

no stress analysis will be performed 
for this element; if MON=l, stresses 
will be calculated. 

Card 2 is repeated NE times. 



Card 2 
Columns 


1-5 

N 

6-15 

X(N,1) 

16-25 

X(N,2) 

26-35 

X(N,3) 

36-45 

X(N,4) 

46-55 

X(N,5) 


Displacement Specification 15, 5F10.0 

Contents 

Node at which displacements are 
specified 

Displacement in direction 1 for node N 

Displacement in direction 2 for node N 

Displacement in direction 3 for node N 

Displacement in direction 4 for node N 

Displacement in direction 5 for node N 


Card 2 is repeated NlD times. 

Cards 1 and 2 are repeated for initial velocity input. 
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IX . IMPACT PARAMETERS 


Card 1 
Columns 
1-5 


6-15 

16-25 

26-35 

36-45 

46-55 

56-65 

66-75 

76-85 


N 

MASS 

K 

N 

M 

I 

C 

Xo 

Vo 


Impactor & Contact Spring Input 15, 8F10.0 

Contents 

The node at which the contact spring 
is connected 

The mass of the impactor 
The spring rate for loading 
The spring exponent for loading 
The spring exponent for unloading 
The permanent deformation exponent 
The permanent deformation multiplier 
The initial position of the impactor 
The initial velocity of the impactor 
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Last Time Step? 


Yes 

Save FILE08 for f 

Restart ^ 
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Figure F.l. CLIP II Flow Chart 
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